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ABSTRACT 

In physically inviscid fluid dynamics, "shock capturing" methods adopt either an arti- 
ficial viscosity contribution or an appropriate Riemann solver algorithm. These tech- 
niques are necessary to solve the strictly hyperbolic Euler equations if flow disconti- 
nuities (the Riemann problem) are to be solved. A necessary dissipation is normally 
used in such cases. An explicit artificial viscosity contribution is normally adopted to 
smooth out spurious heating and to treat transport phenomena. Such a treatment of in- 
viscid flows is also widely adopted in the Smooth Particle Hydrodynamics (SPH) finite 
volume free Lagrangian scheme. In other cases, the intrinsic dissipation of Godunov- 
type methods is implicitly useful. Instead " shock tracking" methods normally use the 
Rankine-Hugoniot jump conditions to solve such problems. A simple, effective solu- 
tion of the Riemann problem in inviscid ideal gases is here proposed, based on an 
empirical reformulation of the equation of state (EoS) in the Euler equations in fluid 
dynamics, whose limit for a motionless gas coincides with the classical EoS of ideal 
gases. The application of such an effective solution to the Riemann problem excludes 
any dependence, in the transport phenomena, on particle smoothing resolution length 
h in non viscous SPH flows. Results on ID shock tube tests, as well as examples of 
application for 2D shear flows are here shown. As an astrophysical application, a much 
better identification of spiral structures in accretion discs in a close binary (CB), as a 
result of this reformulation is also shown here. 

Key words: accretion, accretion discs - equation of state - hydrodynamics: shocks - 
methods: numerical, N-body simulations - stars: binaries: close, cataclysmic variables, 
dwarf novae. 



1 INTRODUCTION 

In both Lagrangian and Eulerian inviscid fluid dynamics, a 
dissipation is normally adopted to handle discontinuities in 
the flow (the Riemann problem). An artiflcial viscosity is 
introduced in SPH, as a shock capturing method, to pre- 
vent particle interpenetration and to smooth out spurious 
heating in the flow to solve the strictly hyperbolic system 
of Euler equations. The introduction of such a small dis- 
sipation is also currently adopted to produce both mass 
and angular momentum transport in SPH phy sically inviscid 
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1997. ll998l ; IOkazaki et al.ll2002l '). Eff'orts were accomplished 
in SPH to solve both the "approximate" and the "exact" 
Riemann problem, either v ia an explicit reformulation of 
the artificial viscosity term (|Balsaralll995l ; lMonaghan|[l997l ; 
iMorris fc Monaghanl Il997 ) or v ia sophisticated Godunov 
algor it hms llYukawa et aL Il997l; [Parshikov 1999; Inutsukp 
2002; 'Parshikov fc Mcdin1 12002 ; 7cha fc WhitworthI |200i 
Molteni fc Bilello..2003 ) which also include an " intrinsic" im- 
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plicit dissipation. In the flrst reformulation of the 

artiflcial viscosity could be necessary because, for "weak 
shocks" or low Mach numbers, the fluid becomes "too vis- 
cous" and angular momentum and vorticity could be non- 
physically transferred. An alternative physical way to solve 
the Riemann problem, based on a reformulation of the EoS 
in the Euler equations, is here presented, where particle SPH 
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pressure terms are recalculated without any artificial viscos- 
ity contribution. 

In this paper we do not solve the Riemann problem by 
proposing a numerical algorithm, i.e. a Riemann solver. In- 
stead we try to solve the Riemann problem, invalidating the 
stability of solutions of the hyperbolic Euler equations of 
inviscid flows, in a strictly physical sense, searching for an 
EoS for ideal gases including the correct dissipative contri- 
bution. This means that, even though we work in SPH for- 
mulation, making the most of our experience, final results 
involve the EoS in general. Since shock flows are non equilib- 
rium events, we pay attention that the EoS: p = (7 — l)pe for 
ideal flows cannot exactly be applied to solve the Riemann 
problem. In fact such an equation is strictly valid only for 
equilibrium or for quasi-equilibrium thermodynamic states, 
whenever the velocity of propagation of perturbations equals 
the sound velocity. Successful results, based on some form of 
mathematical dissipation introduced within the strictly hy- 
perbolic system of non linear Euler equations for ideal flows, 
are obtained due to the necessity to correct such a congeni- 
tal defect. In §2 of this paper, we briefly recall how the SPH 
method works for ideal non viscous flows. In the same §2 we 
also outline which evolution has been accomplished in the 
explicit artiflcial viscosity dissipation description. Instead, 
in §3, we show how to reformulate the EoS, according to the 
Riemann problem for inviscid, ideal gases. Analytical for- 
mulations of how a physical dissipation in non viscous flows 
for irreversible physical processes plays on thermodynamic 
properties is shown in §4, according to statistical thermody- 
namics. In §5 a reflnement of the same EoS is also formulated 
to handle a dissipation that is also right for shear flows. Ap- 
plications to ID shock tubes (|Sod 1978) and to blast waves, 
to solve shocks, to 2D turbulence and 2D shear flows, are 
also shown in §6, as well as a comparison to models adopting 
different dissipation that is also given for the coming out of 
spiral patterns in accretion discs in a CB (see App. A). 



2 THE EULER EQUATIONS AND THEIR SPH 
FORMULATION 

2.1 SPH and artiflcial viscosity for non viscous 
flows 

In the Lagrangian ideal non viscous gas hydrodynamics, the 
relevant equations (Euler equations) are: 



dp 
H 
dv 
Ht 

de 

dt p 
p = (7 - l)pe 
dr 

— — V 

dt 



+ pV v = 

= +* 

P 



continuity equation(l) 
momentum equation(2) 

energy equation(3) 
perfect gas equation(4) 
kinematic equation(5) 



The majority of the adopted symbols have the usual 
meaning: d/dt stands for the Lagrangian derivative, p is the 
gas density, e is the thermal energy per unit mass, $ is 
the external effective force fleld, also taking into account of 
centrifugal and Coriolis contributions. The adiabatic index 



7 has the meaning of a numerical parameter whose value lies 
in the range between 1 and 5/3, in principle. 

The SPH method is a Lagrangian scheme that dis- 
cretizes the fluid into moving interacting and interpolating 
domains called "particles". All particles move according to 
pressure and body forces. The method makes use of a Kernel 
W useful to interpolate a physical quantit y A{r) related to 
a gas particle at position r according to (|Monaghanlll985l . 
il992 ): 



A{r)= I A{r')W{r,r',h)dr' 

' D 



(6) 



W{r,r' ,h), the interpolation Kernel, is a continuous 
function - or two connecting continuous functions whose 
derivatives are continuous even at the connecting point - 
defined in the spatial range 2h, whose limit for /i — !■ is the 
Dirac delta distribution function. All physical quantities are 
described as extensive properties smoothly distributed in 
space and computed by interpolation at r. In SPH terms we 
write: 



A, 



spAj_ 



^A. 



(7) 



where the sum is extended to all particles included 
within the interpolation domain D, rij = pj /rrij is the num- 
ber density relative to the jth particle. W{ri,rj,h) ^ 1 is 
the adopted interpolation Kernel whose value is determined 
by the relative distance between particles i and j. Typically, 
such cubic spline Kernels (|rij |, /i) are in the form: 



1 



ff < \ri 
otherwise, 



s; 1 
«; 2 



(8) 



even though also Gaussian forms are adopted. In ex- 
pression (8), \rij\ — \ri — Vjl/h represents the module of 
the radial distance between particles i and j in units of h. 
Other SPH Kernel analyt ical formulations could be consid- 
ered iFulk fc Quinnlll996l '). However this is beyond the aim 
of this paper because the physical solution of the Riemann 
problem to solve the strictly hyperbolic Euler equations of 
fluid dynamics does not rely on the Kernel choice related to 
the numerical scheme. The Kern el formulation (8) is anyway 
the most popular since the 90's (|Monaghanlll992i 'l. 

In SPH formalism, equations (2) and (3) take the form, 
respectively: 



dvi 
~dt 



+ 



(9) 



(10) 



where Vij — Vi — Vj and rrij is the mass of jth particle. 
For a better energy conservation, the total energy E — 
(e + ^v'^) can also be introduced in the SPH formulation: 

jv ^ \ 
±E, = -J2m,(P^ + P^yS/,W.,+^.-v, (11) 

j=i \ ' ^ / 

of the energy equation: 
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(12) 



In this scheme the continuity equation takes the form: 



(13) 



or, as we adopt, it can be written as: 



(14) 



3=1 



which identifies the natural space interpolation of par- 
ticle densities according to equation (7). 

The pressure term also includes t h e arti fi cial v iscos- 
ity contribution given by Monagha 3 (| 19851 . Il992l ) and 
iMonaghan fc Lattanzid ()l985h . with an appropriate thermal 
diffusion term which reduces shock fluctuations. It is given 
by: 



where 



if Vij ■ Tij < 



otherwise 



(15) 



(16) 



with Cai being the sound speed of the ith particle, 
fij = Vi — rj , <^ , a ~ 1 and / 3 ~ 2. These a and /3 
parameters of the order of the unity (|Lattanzio et al.lll985l ) 
are usually adopted to damp oscillations past high Mach 
number shock front s developed by non linear instabilities 
l|Boris fc Booklll97:j '). SPH method, fike other finite differ- 
ence schemes, is far from the continuum limit. The linear a 
term is based on the viscosity of a gas. The quadratic (/3, Von 
Neumann-Richtmyer-like) artificial viscosity term is neces- 
sary to handle strong shocks. In the physically inviscid SPH 
gas dynamics, angular momentum transport is mainly due 
to the artificial viscosity included in the pressure terms as: 



Ei 



'^^ + 



(17) 



where terms into parentheses refer to intrinsic gas prop- 
erties. 

In SPH conversion (eqs. 9, 10, 11, 13, 14) of mathemat- 
ical equations (eqs. 1 to eq. 3) there are two principles em- 
bedded. Each SPH particle is an extended, spherically sym- 
metric domain where any physical quantity / has a density 



2.2 The Riemann problem in the SPH 
formulation 

Strong shock require a = 1. However for weak shocks and for 
small Mach number fiows, the fiuid becomes "too viscous" 
and both angular momentum and vorticity are transferred 
unphysically. To solve this problem, several solutions are 

proposed: 

a) the formulation of a "limiter" l|Balsaralll995h . multi- 
plying the artificial viscosity terms rji for fij — 0.5(/i + fj), 
where 



(18) 



V ■ i>i| IV X Vi\ + W-*Csi/h' 

It reduces the unphysical spread of angular momentum 
in whirling fiows up to 20 times. It is « 1 for planar shocks, 
while it increases if the tangential kinematics is relevant; 

b) a switch to reduce artificial viscosity 
(|Morris fc Monaghan|[l997l ). In this hypothesis, for each ith 
particle the a evolves according to a decay equation: 



dai 
IF 



(19) 



where a* ~ 0.1, n — h/cs, Cs is the sound speed, and 
the source term Si — fi max( — V ■ Vi,0). 

c) a reformulation of the artificial vis cosity according to 
the Riemann problem (!Monaghanl ll997l ). Particles i and j 
are considered as the equivalent left "F and right "r" states 
of a given contact interface. The ID Riemann problem is 
taken into account along the line joining them. Being the 
Euler equations in conservative form: 



ds df „ 
1 — — — 0, 

dt dx ' 

the simplest Euler technique of integration is: 

s"+^ = s" - [/*(s,,Si_i) - /*(Si_i,Si)] , 

where numerical fluxes (iMartf et al.lll99l] ) 



(20) 



(21) 



(22) 



where the ej are the eigenvectors of the Jacobian matrix 
A = df /ds and A* is an average of A for the "/" and "r" 
states. Alj* are the "jumps" of s across the characteristics: 



Sr — Sl 



2^ Aoj'ej. 



(23) 



For ID ideal flows, the eigenvalues are v, v + Cs and 
v~Cs, where the two including the sound velocity are the ve- 
locities of propagation of sound referred to the frame whose 
fluid velocity is v. Assuming that the jump in the velocity 



profile fW{ri,rj,h) = fW{\ri - rj\,h). Besides, the fluid across characteristics could physically be Vij ■ rij/\rij\ and 



quantity / at the position of each SPH particle could be in- 
terpreted by flltering the particle data for /(r) with a single 
windowing function whose width is h. So doing, fluid data 
are considered isotropically smoothed all around each parti- 
cle along a length scale h. Therefore, according to such two 
concepts, the SPH value of the physical quantity / is both 
the superposition of the overlapping extended profiles of all 
particles and the overlapping of the closest smooth density 
profiles of /. 



that a signal velocity Vgig corresponds to the above eigenval- 



ues, |A*|Aa;* corresponds to Veig 



So doing. 



the artificial pressure contribution in the momentum equa- 
tion is: 



Pi] 



(24) 



where K is an arbitrary constant ~ 1 and pij — 0.5(pi-|- 
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As far as the energy equation is concerned, tiie SPH 
formulation in this case is: 



dt 



N 

i=i 
Kv. 



PiVj 
sig.ijS-ijfij 



+ 



(25) 



being e*j = e* - e*, where e* = O.^ivi ■ rij/\rij\Y + ^i- 
The signal velocities for the ID L agrangian Rieman n 
problem on ide al flows are reported in ( Whitehurstlll99a) . 
on the basis of iGottlieb fc GrothI j 19881 ) and lTord (| 19921 ') 
results. The pressure on the contact interface p* links the 
left and the right states. For systems with one shock and 
one rarefaction wave p* G [pi,pr\ and 



Csl + Csr + {Vl - Vr)i'y - l)/2 



Csl/p\ 



(7-l)/27 



+ Csr/p'. 



(7-l)/27 



27/(7-1) 



(26) 



I n the case of two shoc ks, a more complicated rela- 
tion l|Gottlieb fc GrothI 1 19881 ) also involves the velocity v* 
on the contact interface. However, for some practical pur- 
poses, p* = max(p;,pr) and v* = 0.5{vi + Vr). Notice that 
in the presence of two very strong rarefaction waves p* — 0. 

The two Lagrangian velocities of waves spawned by the 
interface are: 



Vl 



A, 



Vl - Csl 



and 



Vr + Csr\l + 



Vr + Cs 



27 

if P* /pi 1 

(7-i)(pVp.-i) ] ifp7p^>i 



27 



(27) 



(28) 



if p* /pr 5? 1 



In the Lagrangian description, being Vsig the speed of 
propagation of perturbation from i to j and vice-versa {I -H- 
r), 



— Csj — Vj • 



Csi + Csj — Vij ■ f^ij/\f^ij\ 



(29) 
(30) 
(31) 
(32) 



having considered the versus goi ng from i to j. By per- 
forming some numerical experiments. iMonaghanI ([1993) also 
suggested other similar algebraic expressions. 



2.3 Some cautions in using artificial viscosities 

In[sh^ l|l992h . some cautionary remarks are reported on the 
adoption of artificial dissipation. In particular, "Because the 
magnitude of the viscous term does not affect the net shock 
jump conditions, many numerical schemes implicitly or ex- 
plicitly incorporate the trick of artificial viscosity for halting 
the ever-growing steepening tendency produced by nonlin- 
ear effects, thereby gaining the automatic insertion of shock 
waves wherever they are needed (in time-dependent calcu- 
lations)." However, the numerical viscous term should be 
large enough to spread out shock transitions only over a few 



resolution lengths, making shock waves resolvable. In this 
sense, any mathematical dissipation should be considered a 
useful mathematical "trick" to get correct results. 



3 HOW EOS MATCHES THE RIEMANN 
PROBLEM 

The solution of the Riemann problem is obtained at 
the interparticle contact points among particles, where 
a pressure and a velocity, relative to the flow dis- 

co ntinuity, a r e computed. T his is also clearly shown 

in 'Parshi kovl (Il999l): llnu tsuka' ("2002'); 'Parshikov fc Medin| 
(2002); C ha fc Whitwort h (2003); Moltcni fc BilcUo (2003j), 
where the new pressure p* and velocity v* are reintroduced 
in the Euler equations instead of p and v to obtain the 
new solutions compatible with inviscid flo w discontinuities . 
In SP H, wc pay attention in particular fParshikov 'l999|; 
Inutsu ka 2002; Parshikov & Mcdin 2002; Cha fc Whitworty 
20031 ) to the particle pressures Pi and pj , in the SPH formu- 



lation of the momentum and energy equations (9) and (10 
or 11), whose substitution with pressures p* and p*, solu- 
tions of the Riemann problem, excludes any artificial viscos- 
ity adoption, although a dissipation, implicitly introduced, 
is necessary. Therefore, the solution of the Lagrangian Rie- 
mann problem, only in its pressure computation, is enough 
to interface SPH with the Riemann problem solution. 

The key concept is that the physical interpretation of 
the application of the artificial viscosity contribution in the 
pressure terms corresponds to a reformulation of the EoS for 
inviscid ideal gases, whose equation: 



P\gas = (7 - l)pe 



(33) 



is strictly applied in fluid dynamics when the gas com- 
ponents do not collide with each other. In the case of gas 
collisions, it modifles in: 



p* = {-f — l)/9e -I- other. 



(34) 



The further term tak es into account the velocity of 
perturbation propagation (|Monaghanlll997l ). This velocity 
equals the ideal gas sound velocity Cs whenever we treat 
non coUisional gases in equilibrium or in the case of rarefac- 
tion waves. Instead, it includes the "compression velocity": 
Vij ■ Vij /\rij I in the case of shocks (eq. 29). In the first case, 
we write the EoS for inviscid ideal gases as: 

p\aas = -cl, (35) 
7 

where Ca = (tipI pY^'^ ~ Ylk~1 ~ l)^]"'"''^- Instead, in the 
second case, the new formulation for the EoS is obtained 
squaring Vsig,i^j, so that Cs(l - Vshockjcsf is an energy per 
unit mass in the case of shock. Hence: 

\2 



(1 



if Vs 



< 



if V shock > 

In the SPH scheme, being: 



Pi = 



7 V 



T-^shock.i 



if v,i 



otherwise. 



(36) 



(37) 



(38) 
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This formulation introduces the "shock pressure term": 
pi'^lhock ^ 2vshockCs) whose Unear and quadratic power 
dependence on Vij ■ rij/\rij\ is analogue to both the lin- 
ear and the quadratic components of the artificial viscos- 
ity terms (15). The linear term oc CsVshock is based on the 
viscosity of a gas. The quadratic term oc {csVahock)'^ (Von 
Neumann-Richtmyer-like) is necessary to handle strong 
shocks. These contributions involve a dissipative power, 
whose effect corresponds to an increase of the gas pressure. 
Therefore, we adopt the formulation (eq. 21) as p* and p* in 
the SPH formulation of the momentum (eq. 9) and energy 
equations (eqs. 10 or 11): 



dvi 

Hi 

dt 



N 
N 



dt 



N 



(39) 



(40) 



(41) 



This simple reformulation, allows us to keep the same 
Courant-Friedrichs-Lewy condition as for the timestep com- 
putation, substituting only the sound velocity Cs with Cs — 

^ shock . 

Notice that the comparison of physically dissipative 
shock pressure terms pivg^ock ~ '2vshockCs)/y with those rel- 
ative to SPH artificial viscosity: Pilgasilij gives for a and /3 
parameters the physical equivalence for the ith particle: 



P = 



^sij I'^ij I 



(42) 
(43) 



where Csij = 0.5(csi + Csj). 

The linear and the quadratic dependence of a and /? 
on Csij\rij\/{csih) perfectly correlates the physical dissipa- 
tion of a perfect gas both to a bulk viscosity and to a 
quadratic Von Neumann-Richtmyer-like dissipation, able to 
handle strong shocks, avoiding any dependence on h, con- 
trary to SPH artificial viscosity. In fact, both h and h'^ alge- 
brically erase applying eqs. (15) and (16). Compared to SPH 
constant a ~ 1 and /3 ~ 2, such correlations show that a and 
P equal zero for \rij \ = and that a = 2 and /3 = 1 in a ho- 
mogeneous and isotropic gas for particle mutual separations 
\rij\ = h. Both parameters linearly and quadratically vari- 
ate, respectively, varying either \rij\/h, or Csij/csi, or both. 
Hence, both a and /? decrease for closer particle neighbours 
and/or for colder companions, whose Csj is smaller than Csi- 



4 CONSEQUENCES OF A SHOCK PHYSICAL 
DISSIPATION ON THE 

MAXWELL-BOLTZMANN STATISTICAL 
THERMODYNAMICS 

Thermodynamic properties of a system in ther mal equilib- 
rium are describe d by the partition function ifplleif f 1965)1 : 
iMcClelland Ill973l ) Z: 



energy levels 



quantum states 



E 



complexions 



E 



-l3Uj 



(44) 



The most probable distribution of systems in the ensem- 
ble with energies Ui and with quantum states of the system, 
given by the Maxwell-Boltzmann law, are: 

— ^ -a -BUi 



and 



(45) 



(46) 



respectively, where e~" = Z/N, being N = . rii the 
total number of systems in the ensemble, U = UiUi the 
total energy of the ensemble and being d the degeneracy 
of the energy level Ui. The parameter /? must be strictly 
positive to have Z meaningful. /3 = {KbT)~^, where Kb 
is the Boltzmann constant and T is the temperature. Being 
Ui classically proportional to the square of the linear mo- 
mentum: {Ui oc v'f) for non interacting free atoms, without 
considering their internal energy levels, each exponential in 
the summations is related to a Gaussian distribution func- 
tion. 

Since in fiuid dynamics equilibrium conditions are a spe- 
cial case, the application of thermodynamic laws is an ev- 
ident approximation, in particular whenever non reversible 
events occur like shock waves. We assume the hypothesis 
that even in conditions of a fluid motion, the Maxwell- 
Boltzmann law holds, where the net effect of any physi- 
cal dissipation is to reduce the peak value of the Maxwell- 
Boltzmann distribution, enlarging its spread, conserving 
both N and U . Hence, any physical dissipation is introduced 
as a < -D ^5 1 term [D = 1 involves no dissipation) mul- 
tiplying pUi or pUj in the above exponentials as: fiDUi or 
PDUj, so that: 



energy levels 



quantum states 



Z = 



E 



-flDUi 



E 



complexions 



E 



(47) 



This is necessary to determine the desired result, with- 
out any alteration of the system's degeneracy, which can- 
not be considered in so far as we are dealing with an ideal 
gas. In a microscopic description, physical dissipation irre- 
versibly converts ordered kinetic energy into chaotic kinetic 
energy. That is the conversion of macroscopic kinetic en- 
ergy into thermal energy, in so far as the thermodynamic 
system is always in thermal equilibrium, obeying Marko- 
vian statistics. Atoms or molecules do not interact with 
each other, and their internal quantum energies do not af- 
fect the global thermodynamics. In the Maxwell-Boltzmann 
distribution hypothesis, this implies a transition from one 
Maxwell-Boltzmann statistical distribution to another one, 
whose half weight at half height is larger. If such a transi- 
tion occurs during the shock crossing, the whole thermody- 
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namic system could make an instantaneous transition be- 
tween two thermodynamic equilibrium states. Hence, the 
larger the half weight at half height of the statistical distri- 
bution, the larger the physical dissipation occurred. Thus, 
writing 13* = D{KBT*y^ or, that is the same: /3* = D/S = 
D{KbT)-\ 



physically corresponding to the EoS for ideal flows, even 
taking into account a shock occurrence, exactly in the form 
(36, or 37-1-38) previously written. Enthalpy and Gibbs free 
energy are given by: H* = U* + p*V and G = U* -\- p*V — 
TS* , respectively. 



energy levels 



complexions 



quantum states 



3 



(48) 



The net effect of such a dissipation is the desired en- 
largement of the width at half height of each Gaussian dis- 
tribution for Z, conservi ng both N and U. Thermodyna mic 
properties are obtained jRl leif f 1965)1 : lMcClellan(nil973l ') in 
the same way as they are currently given using /3* instead 
of /3 as: 



• Internal Energy: 
din Z/dT\v 



U* = 



KgT 



KbT 



dP*/dT\v \ D l3dD/dT\v 
keeping fixed the volume V; 

• Entropy: 

From the Boltzmann law, 

KbN U* 

= — m Z H ; 

D T ' 

• Free Helmholtz Energy: 
NKbT 



dlnZ 



dT 



(49) 



F=U*- TS* = 
• Pressure: 



D 



InZ: 



(50) 



(51) 



In the hypothesis that dissipation does not explicitly 
depend on V, 



P = 



dF* 
~dV 



N K bT In Z 



D 



V 



(52) 



In the classical case of free, non interacting atoms, 
din Z/dV\T = 1/V, hence 



P = 



p* ^ N 



iV j_ 

Thus, being /3* = = D{KbT)-\ 
KbT 1 



V D' 



(53) 



(54) 



Assuming D = 1/(1 - v shock, i/csiY ^ 1; Vshock,i = 
'"i.j ■ fij/lfi.jl if ""iJ ■ '''i.j < 0, Vshock,i = Otherwise; we 
write an expression that does not depend on volume V and 
corresponding to a function that, even depending on T, it 
can be easily handled by considering the derivatives in (54). 
Thus, the EoS is: 



V 



-Y if Vshock < 0, i.e. D <l 



(55) 



if Vshock ^ 0, i.e. D = 1 



5 WHICH EOS IS ALSO RIGHT FOR SHEAR 
FLOWS? 

To make a generalization, we need only one general EoS 
and not a separation of the EoS according to the kinematic 
of the flow. To this purpose, we can generalize the EoS: 
p* = pcs(l - Vshock/csY h as: 



P 2 /, r^-VR 
p = -c^ 1 - C — 
7 V Cs 



1 

where C 



1 for vr 



(56) 

j/\rij \ < 0, whilst C -> 



otherwise. A simple empirical formulation can be. 



(57) 



C = — arccot 

TT 

where i? 3> 1. -R is a large number describing how much 
the flow description corresponds to that of an ideal gas. To 
this purpose, R « X/d, being A oc p~^^^ the molecular mean 
free path, and d the mean linear dimension of gas molecules. 

Although the physical meaning of vr in the field of a 
free Lagrangian particle technique is clear, it could be con- 
troversial in an Eulerian description. In this case, the local 
physical properties of medium have to be taken into account 
because it is necessary to correlate the approaching of two 
particles to a local compression. To convert expression (56) 
to a more general form, we pay attention to the continuity 
equation regarding the numerical concentration n as: 

dn 

— — h n V ■ V ■ 

dt 







(58) 



to determine vr as: 



VR = 



-4/3 



dn 



dt 3^^ dt 3" 



-1/3. 



so that eq. (56) can also be written as: 



7 

where 



-1/3 



C- 



V • V 



3cs 



C=iarccotf7?!i::^ 

TT V 3Cs 



(59) 



(60) 



(61) 



EoS in the form of eq. (36, or 37-1-38, or 56, or 60) are 
equivalent and effective strictly to solve the Riemann prob- 
lem. However, if there is a velocity shear in the inviscid ideal 
flow, the physical dissipation within the EoS (eq. 36, 37-1-38, 
or 56) is activated, despite the lack of any shock, whenever 
close portions of fluid approach each other, showing different 
densities. This is obviously due to the fact that the physical 
dissipation in the EoS is activated in the approaching of only 
two local particles with each other. Looking at the correla- 
tion of physical dissipation with a and /? SPH parameters 
(eqs. 42, 43), a dissipation still occurs even though dissi- 
pation depends on local thermodynamic conditions. This is 
an initial advantage compared to classical SPH. Such a dif- 
ficulty, more limited compared to that arising in SPH (or 
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whatever is the technique adopting an explicit artificial dis- 
sipation), is totally absent in the form of eq. (60), because 
it is free of any dissipation whenever V ■ i) ^ 0. As an exam- 
ple, any simulation relative to non viscous Keplerian fiows, 
where V • i; = at the beginning, does not produce any dis- 
sipation and fiuid motion is stationary endlessly. This is a 
necessary introduction to the tests shown in the subsequent 
sections. In its form (60) the reformulated EoS for non vis- 
cous ideal gases includes a physical dissipation, depending 
on general local conditions, which is activated only when a 
local compression occurs. Hence, dissipation in SPH-based 
techniques is no longer due to the approaching of only two 
particles, but to the general physical properties of the local 
medium, where a lot of particles have to be considered. 

According to such second reformulation of the EoS, also 
right for shear flows, the algebraic relations determining the 
a and /3 parameters (similarly as for eqs. 42 and 43) are: 



= 2- 



-1/3 VT 

V ■ Vj 

h Vij ■ Ti 

-2/3 



(62) 



(63) 



Such relations compare with eqs. (42, 43), respectively, 
in the case of a pure R iemann problem in fiuid dynamics 
ijSieder fc Riffertll2000l l . being V ■ ~ v^J ■ rijh/\rij\^. 



6 TESTS 

In this section we compare solutions in an SPH approach, 
where the modified EoS is taken into account without any 
thermal conduction contribution, with analytical solutions, 
whenever it is possible, as well as with those in SPH where, 
typically an expli cit thermal conductive term is also included 
ijMonagha ^ Il992l ). Tests regard typical flow cases where ei- 
ther collisions (shocks), or turbulence, or shear flows are 
involved. Throughout the simulations, the initial particle 
smoothing resolution length is ft = 5 ■ 10~ whilst the di- 
mension of the whole computational domain D is chosen to 
detect particle turbulence, if it exists, being D/h ^ 100 in 
the most of discussed cases. The adopted 7 value is 7 = 5/3, 
whenever not explicitly specifled. Results concerning the 
SPH4-E0S approach in the form of eq. (60) will be clearly 
discussed whenever substantial differences develop. 



6.1 ID Sod shock tubes 

The behaviour of shock waves moving in the prevailing 
flow is analyticall y described by the Rankine-Hugoniot 



"jump conditio ns" l|Leveauelll99a lHirsch|[l997l : lTordTl999l : 
lBatchelo? '2000'') . These are obtained by spatially integrating 
the ID hyperbolic Euler equations across the discontinuity 
between the two flow regimes left-right in their Eulerian for- 



(64) 
(65) 
(66) 



malism: 




dp 


d 


dt 


dx 


dpv 


d 


'df 


dx 


dpE 


d 


at 


dx 



where E = v'^ + whose conservative analytical form 
can be synthesized as: 

In the limit of zero thickness of the shock discontinuity, 

S{WI - Wr) = f{wi) - f{Wr). (68) 



Under these conditions a requirement for a unique 
single- valued solution is that the solution should satisfy the 
Lax entropy condition (|Leveauel Il992l : [Hirschl Il997l : iTord 
119991 ) f'(wi) < s < f'{wr), where f'{wi) and f'{u)r) are 
the characteristic speeds at upstream and downstream con- 
ditions, respectively. The integrated form of the Rankine- 
Hugoniot jump conditions are: 



Pm — PrVr 

..2\ I ..2\ 



S{pi - pr) 
s{plVl - PrVr) = (Pivf) - (pr^r) 
s{piEl - prEr) = pmEl - prVrEr, 



(69) 
(70) 
(71) 



It is shown, after some algebraic passages, that the 
shock speed is: 



S = Vl+ Csl 



i + l+l El^i 



27 



Pi 



1/2 



(72) 



where Csi = {iPi / piY^"^ ■ In the case of stationary shocks 
being both the upstream and downstream pressures positive, 
there is an upper limit on the density ratio: pi j pr ^ (7 + 
l)/(7 — 1). However, this limit is currently applied also to 
non steady shocks. 

In this section a comparison of analytical an d our SPH - 
Riemann (SPH-l-EoS) ID shock tube test resuhs (|Sodlll978l ). 
also with the initial particle configuration (time T = 0), is 
made. Notice that the so called analytical solution of the 
Riemann problem is obtained through iterative procedures 
left-right to the discontinuity using the Rankine-Hugoniot 
"jump" solutions. Figg. 1 and 2 show results concerning the 
particle density, thermal energy per unit mass and velocity, 
after a considerable time evolution at time T = 100. The 
whole computational domain is built up with 2001 particles 
from X = to X = 100, whose mass is different, according 
to the shock initial position. At time T = all particles 
are motionless, while the ratios pi/p2 = 3 and ei/e2 = 2 
in Fig. 1, and pi/p2 = 3 and ei/e2 = 1 in Fig. 2, between 
the two sides left-right. The first 5 and the last 5 particles 
of the ID computational domain, keep fixed positions and 
do not move. The choice of the final computational time is 
totally arbitrary, since the shock progresses in time, v = 
at the beginning of each simulation. Hence the adimensional 
temporal unity is chosen so that J^dx/cs = 1. Being the 
sound velocity initially constant, this mathematically means 
that Jy dx = Cs, so that I — Cs- 

Results, where our SPH-Riemann (SPH-fEoS) solution 
is applied to SPH, display a good comparison with the 
analytical solution. Discrepancies involve only 4 particle 
smoothing resolution lengths at most, except for numeri- 
cal solutions corresponding to analytical vertical profiles re- 
garding thermal energy where, as SPH solutions, discrep- 
ancies are larger. This means that, according to the cau- 
tionary remarks in §2.3, the physical dissipation introduced 
in the EoS (eqs. 36, or 37-1-38, or 56, or 60) is effective. 



8 G. Lanzafame 



0.03 



0.02 



> 0.01 


-o.os 

0.04 
0.03 
0.02 
0.01 

0.0040 
0.0035 

0.003 
0.0025 

0.002 
0.00120 

15 

10 



^ T = 100 - SPH 



' I ' ' ' ' r 

T = 100 - SPH 



^ T = 100 - SPH 



T = 100 - SPH 



-10 




X 



H \ h 



- T = 100 - SPH + EoS 



n I I r 



n \ I r 



H — \ — \ — h 



T = 100 - SPH + EoS 

H — \ — \ — \ — \ — \ — k 



T = 100 - SPH + EoS 

H — \ — \ — \ — \ — \ — k 



T = 100 - SPH + EoS 



J I I L 



10 



-10 




X 



H — — \ — h 



H — — \ — h 



H — — \ — h 



J I I L 



10 



Figure 1. ID shock tube tests as far as both analytical (solid line) and our SPH-Riemann (SPH+EoS, dots) results are concerned (right 
side plots). Density p, thermal energy e, pressure p and velocity v are plotted at time T = 100. Density and thermal energy of particles 
initially at rest at time T = refer to values plotted at the two edges of each plot. The initial velocity is zero throughout. SPH results 
are also reported (left side plots). 



Any arbitrary artificial viscosity contribution is totally ab- 
sent, as well as any thermodynamic properties of neigh- 
bour particle are take n into account, as t he applicat ion of 
Godunov-type solvers (|Yukawa et al.l 1997;: Parshik^ 19991: 
Inutsuka"2002': 'Parshikov fc Medinll2002l : ICha fc WhitworthI 
2003 : Molteni & BilcUo 200^) does. InSra+EoS approach, 
a less evident Gibbs phenomenon (|Gibbs|[l898al lbl). up to 
20 -T- 30% less affects the numerical solution. 

To check the reliability of the adopted EoS (eqs. 36, or 
37-1-38, or 56, or 60), the same ID shock tube tests were 



performed adopting a smaller particle smoothing resolution 
length h = 5 - W~^, working with 20001 particles and adopt- 
ing the same initial and boundary conditions. This check is 
fundamental to prove that particles do not interpenetrate 
with each other and to verify the final results as a function 
of the improved adopted spatial smoothing resolution. Re- 
sults are shown in Figg. 3 and 4, where the comparison with 
both previous results and with the analytical solution are 
also reported. If results are good enough in a low-medium 
smoothing resolution, they are even better in higher spa- 
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Figure 2. The same as in Fig. 1. In this example, the initial discontinuity does not ailect the thermal energy per unit mass, being the 
gas initially isothermal, while the initial velocity is zero throughout. 



tial resolution, confirming that the shock profile depends on 
spatial resolution length. However this is a natural result 
emerging whichever is the adopted numerical technique. 

For the SPH+EoS modelling where eqs. (36, or 37+38, 
or 56) are used, the SPH a and 13 parameters of artificial 
viscosity counterpart (eqs. 42, 43) are shown in Figg. 5 and 
6 for both Sod shock tube tests. These calculations must be 
considered for each neighbour of each particle because they 
take into account both Csi/csij and \r\ij/h, as well as both 
(?si/<?sii and |r|fj//i^. In the fiat motionless zones, far from 
the shock discontinuities, \r\ij = h, so that a = 2 and /? = 1, 



whilst those neighbours whose \r\ij = 2h are not taken into 
account because SPH interpolations are ineffective at the 
Kernel edge. In these motionless regions, despite a = 2 and 
P = 1, any physical dissipation in SPH+EoS is not active, 
as in SPH, because Vij ■ Vij = 0. Notice that, according to 
the EoS eqs. (36, or 37+38, or 56), whenever nj ■ Vij ^ 0, 
any dissipation is physically prevented. Thus, whichever are 
the a and /3 counterparts, their role is useless. Both a and 
P assume multiple values wherever any kinetic evolution of 
the flow is recorded, especially in the proximity of the dis- 
continuities, not only because of variations of particle mu- 
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Figure 3. Enlargement of discontinuities of Fig. 1 botii for 
h = 5 ■ 10-2 (long dashes) and for ft = 5 • IQ-^ (dots) SPH- 
Riemann (SPH+EoS) resolution lengths. The analytical solution 
is also plotted (solid line). 
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Figure 4. Enlargement of discontinuities of Fig. 2 both for 
h = 5 ■ 10-2 (long dashes) and for ft = 5 ■ 10"'^ (dots) SPH- 
Riemann (SPH+EoS) resolution lengths. The analytical solution 
is also plotted (solid line). 



Figure 5. Parameters a and /3 of the artificial viscosity counter- 
part for the SOD ID shocktube test of Fig. 1. 



Figure 6. Parameters a and /3 of the artificial viscosity counter- 
part for the SOD ID shocktube test of Fig. 2. 



values. Besides, as Figg. 5 and 6 show, any smoothing on 
the vertical discontinuities reflects the SPH smoothing of 
vertical linear tracks in the shock profiles. 



tual separation and sound velocity, but also because of small 
statistical fluctuations in the particle kinematic and thermal 
properties. This shovsrs both the tendency for the closer inter- 
acting particles (smaller \r\ij/h < 1) to have smaller a and 
P values - whilst it is the opposite for the farther interacting 
particles (larger 1 < \r\ij/h < 2), in so far as Vij ■ Vij < 
- and very small horizontal "dual" behaviour of a and /3 



6.2 ID Blast wave 

Whenever in a shocktube the ratios pi/pi — ei/e2 3> 1, and 
consequently pi/p-z ~ 1, and vi = V2 = 0, such a discon- 
tinuity is called a "blast wave". Being v = at the begin- 
ning of each simulation, the adimensional temporal unity 
is chosen as previously written in the ID Sod shocktube 
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Figure 7. ID blast wave test as far as both analytical (solid line) ajid our SPH-Riemann (SPH+EoS, dots) results are concerned (right 
side plots). Density p, thermal energy e, pressure p and velocity v are plotted at time T = 5. Density and thermal energy of particles 
initially at rest at time T = refer to values plotted at the two edges of each plot. The initial velocity is zero throughout. SPH results 
are also reported (left side plots). 



test before. Figg. 7 and 8 show a comparison of SPH and 

SPH+EoS results with the so called analytical solution, as 
shortly described in the previous subsection, after a consid- 
erable time evolution at time T = 5. However, especially in 
these cases, the analytical solution is considered corrected 
in so far as pi/p2 ^ (7 + l)/(7 — !)• In the blast wave 
test here considered, pi/pi — €i/e2 ~ 10*, while other spa- 
tial, initial and boundary conditions, as well as the particle 
spatial smoothing resolution length arc identical to those 
chosen in the previous test. Figg. 7 and 8 show that SPH 



and SPH+EoS results globally compare with each other and 

that they also compare with the analytical solution wher- 
ever pi/ p2 ^ (7 + 1)7(7—1), that is wherever the Rankine- 
Hugoniot jump conditions hold. Beyond this limit, even the 
so called analytical solution is considered incorrect. Being 
7 = 5/3, the comparison is meaningful within pi/p2 ^ 4. 
However, the SPH+EoS modelling has the advantage of 
offering a better solution, compared to the SPH one, for 
pi/p2 > 4 (and therefore for pi /p2 > (7+1)7(7— 1)) as Fig. 8 
clearly shows, where the SPH numerical solution suffers from 
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smoothing resolution hides this effect because of the higher 
artificial damping due to a higher particle smoothing resolu- 
tion length h (eq. 15-16). Moreover, in SPH, even the choice 
of the arbitrary parameters a and /3 should be linked to the 
specific physical problem. Instead, in SPH-l-EoS, the damp- 
ing is strictly physical and local and any "blimp" effect is 
strongly reduced, especially for ID blast waves, where strong 
discontinuities in the flow deeply affect the SPH numerical 
solution producing intrinsic numerical instabilities close to 
the propagating discontinuities (Figg. 7-8). The higher the 
spatial smoothing resolution (the smaller h), the higher the 
"blimp" instabilities. 

Fig. 9 shows the a and /3 parameters of the artificial 
viscosity counterpart for the ID blast wave of Fig. 7. Com- 
putational considerations for both parameters are equivalent 
to those written in the previous ID Sod shock tube test. No- 
tice that in such a test, in the denser region of the blast wave, 
the number of neighbour particles involved in the computa- 
tion of such parameters is larger than in the previous test, 
as evidenced by the multiple a and /3 values. 



Figure 8. Enlargement of previous Fig. 7 showing SPH instabili- 
ties at contact discontinuity. The Eos-(-SPH profile (both in dots) 
does not of such instabilities. The analytical profile (solid line) is 
also shown. 



Figure 9. Parameters a and /3 of the artificial viscosity counter- 
part for the ID blast wave test of Fig. 8. 



some instabilities. Moreover, the SPH-l-EoS solution does 
not seem to suffer from any "blimp " effect at the con tact dis- 
continuity, as briefly discussed in iMonaghanI l| 19971 ). where 
a modiflcation of the artificial viscosity term is proposed 
(see also §2.2), as far as the velocity profile is concerned. 
Such effect comes out whenever a spatial high smoothing 
resolution is working together with an explicit handling of 
dissipation through an artificial viscosity damping to solve 
the Riemann problem of flow discontinuities. A low spatial 



6.3 2D radial spread and migration of an annulus 
Keplerian ring 

The 2D radial spre ad and m igration of an annulus ring is 
widely described in Frank et al. (2002) in the case of a con- 
stant physical viscosity i/. At time T = 0, the surface density, 
as a function of the radial distance r, is described by a Dirac 
5 function: E(r, 0) = M5{r ~ ro)/2nro, where M is the mass 
of the entire ring and ro is its initial radius. The following 
time, the surface density is computed via standard methods 
as a function of the modified Bessel function Iin(z): 



Ti{x, t) — -^r ^3 



-1/4, 



-h/i{2x/T) 



(73) 



where x — r/vo and r = 12vT/ro^, expressing ra- 
dial distances in ro units and time in viscous dissipation 
time units, respectively. J^Tj{x,T)dS — 2-k ^ J^{x,T)dr = 
const equals the annulus mass throughout. Time is normal- 
ized so that r = 1 corresponds to the Keplerian period 
relative to the ring at ro = 1. Exa mples of SPH sprea d 
on this argument can be found in iFlebbe et al.l l|l994 ): 
Soeith fc R iffcrt (l999h: ISpeith fc Kiev! (|2003 '). as weU as in 
Costa et al. (2009.) in SPH physically inviscid hydrodynam- 
ics on the basis that the SPH shear dissipation in Eulerian 
non viscous flows can be compared to physical dissipat ion 
( Molt eni et aD Il99ll : iMurravl Il996l : lOkazaki etall |2002| ) in 
a full Navier-Stokes approach. I n particular an exhaust ive 
comparison can also be found in lLanzafam3 (|2008l . l2009l l. 

A signiflcant comparison of SPH-|-EoS (eqs. 36, 37-1-38, 
or 56, not eq. 60) to SPH is shown in Fig. 10, where XY den- 
sity contour map plots are shown at the same times, to show 
whether dissipation in SPH+EoS approach gives results bet- 
ter fitting to the analytical viscous solution {v « Csh) than 
the classical SPH dissipation. The radial distributions of sur- 
face density are shown in Figg. 11-12, as well as the radial 
profile of the theoretical surface density, according to the re- 
stricted hypotheses of the standard mechanism of physical 
di ssipation (constan t dissi p ation, zero initi a l thick ness). As 
in ISpeith fc Rifei^ ()l999l ): ISoeith fc KIctI (|2003t ). the ini- 
tial ring radius is at ro = 1, whose thickness is Ar = 0.5, 
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Figure 10. XY plots of ring density contour maps. Times are 
reported for each configuration both theoretical and numerical 
(SPH or SPH+EoS). 



is composed of 40000 equal mass (m, = 2.5 • 10"^^) pres- 
sureless Keplerian (1; = VKepi, V ■ w = at T = 0) SPH 
particles, with h = 9 ■ 10~^, with Cs — 5 ■ , and with 
initial density radial distribution corresponding to the ana- 
lytical solution at time T, whose r = 0.018. To this pur- 
pose, a random num ber generator has been used, as in 
ISoeith fc Klevl (|2003h . The central accretor mass is normal- 
ized to M = 1. The kine matic s hear dissipation is estimated 
as V PC c,h (iMolteni et a l. 1991: Murravlll996l : IOkazaki et al.1 
l2002l : lLanzafamdl2q08l. 1200a'). SPH results o f this test com- 
pare w ith th ose of ISpeith fc RiffeTtl (|l999l ): ISpeith fc Klevl 
(|2003l ') and of lCosta et al.l (|2009l '). Being in an Eulerian fluid 
dynamics, in SPH the shear effect is produced by the arti- 
ficial viscosity, and the consequent artificial pressure terms, 
having removed the physical pressure terms. ISpeith fc Klevl 



discussed its limits, treating a ring free of any pres- 
sure force, where a viscous stress tensor similarity of vis- 
cosity is considered. In SPH-l-EoS the physical dissipation, 
appearing as further terms in the EoS, determines the shear- 
ing effects, which is still effective also removing the physical 
pressure terms as in the SPH counterpart, ft is decisive the 
fact that the SPH+EoS radial density much better fit both 
the theoretical radial profile and its radial migration towards 
the central accretor. Notice that these results are obtained 
according to the cubic spline (8) Kernel analytica l expres- 
sion and according to the correlation v ^ Csh (Molte ni et al.l 
Il99l[ ) in the expression where r = VluT /r^ ■ 

The arbitrary assumption of a different p would involve 
the appropriate time T to get the same r. In fact, the vari- 
ation of the initial h, and /or Cs, to get a smaller (larger) 
u implies a longer (shorter) time T to get the same tau to 
produce the same radial distribution of particles. 

Instead, different radial distributions are obtained as- 
suming a different interpretation of the analytical expres- 
sion for V. In fact, considering v ~ Q.laspHCsh (|Murravl 



Figure 11. Surface density E(r, t) in 10"^^ units as a function 
of radial distance from initial configurations at r = 0.018 when 
radial profile equals the theoretical analytical one. Subsequent 
snapshots are reported for each configuration both theoretical 
(left side plots) and numerical (SPH or SPH+EoS). 




Figure 12. As Fig. 11 for more evolved times. 



1 19961 : lOkazaki et aLlliool ). with asPH « 1, it is necessary 
to perform the numerical simulations for a time T ten times 
longer to get the same r, in spite of the initial h and Cs axe 
unchanged. This involves a too large annulus spread as far 
as the numerical results are concerned in both cases, com- 
pared to the anal ytical solution. To this purpose, even the 
successful test in llVladdison et al.1 (|l996l ^: iMurravl (jl996h is 
not conclusive because in those papers the comparison with 
the analytical solution is successful up to a time T = 0.8 



14 G. Lanzafame 



which, in that parametrization, would involve a too short 
final T ~ 3.32 • 10"*, that is much shorter than our r scales. 
This is also demonstrated in their results because, together 
with their annulus spread, a very short migration of the en- 
tire annular distribution towards r = is also shown as a 
consequence of the fact that the final ring configuration is 
due to a too short tim e evolution. Instead, S P H results of 
ISpeith fc RifertI (| 19991 ) : ISpeith fc Kiev! jiooj ): ICosta et all 
compare with our SPH results because the spread 
and migration is not meaningless. It is the SPH+EoS that 
is better in so far as v ^ Csh is considered. 

a and j3 artificial viscosity parameters of the SPH coun- 
terpart, for SPH+EoS modelling according to eqs. (42, 43), 
are simply expressed by a = 2rij/h and j3 = (rtj/h)^, since 
the whole annulus ring is isothermal in such shockless and 
pressureless simulations. This means that, in their radial 
profile - here not shown for the sake of simplicity - both 
parameters span from small positive values in the denser 
regions of the annulus ring, up to their maximum values 
(q = /? = 4), at the beginning of each SPH-|-EoS model. 
However, their minimum values increase, towards their max- 
imum limit (4), because of the increase of the particle mean 
free path, as time goes by. 

Circular rings, appearing in Fig. 10 for both numerical 
schemes, are an unavoidable effect due t o the Lagrangian 
part icle-base d techniq ue, as discussed in ISpeith &: RiffertI 
}l999); Spci th fc KlevI ^003) ■ This effect cannot be present 
throughout in the theoretical XY plots because of the 
random-based representation of data points. 

Instead, if the strictly annulus ring has to be strictly 
Keplerian, being \7 ■ v — 0, the SPH+EoS form (60) does 
not produce any radial transport, being zero the intrinsic 
physical dissipation. Hence, the same ring configuration at 
T = 0.018 stays endlessly the same. Hence, results on the 
radial spread and radial migration of such a test are not 
shown. 

6.4 2D slipping flow of a single row of particles 
within a horizontal bounded pipe 

In this test we compare results concerning the ideal flow 
of one horizontal row of particles in a horizontal bounded 
pipe. To this purpose, nine horizontal rows of equal mass 
particles are considered, whose m = 10~^. The computa- 
tion is a fully 2D computation. The four up and the four 
down external rows are densely populated of permanently 
isothermal fixed particles, whose mutual horizonthal sep- 
aration equals /i/4, whose thermal energy per unit mass 
equals that of the flowing particles within the pipe. The in- 
termediate flowing row of particles consists of equally spaced 
particles whose mutual separation equals h and whose uni- 
form velocity v = h = 5 ■ 10~^. The initial Mach number is 
~ 3.5- 10~^ throughout. h/A is also the mutual vertical sepa- 
ration throughout. Our purpose is to determine a permanent 
condition where \7 ■ v — because the analytical solution 
for an ideal gas is that of a constant flow, whose velocity 
and whose thermodynamic properties are constant in time. 
To this purpose, the denser are the boundaries (in terms of 
numerical density), the better are such conditions since V-v 
statistically better equals zero. In such a test, the analytical 
solution is perfectly known. However, our purpose is to verify 
whether SPH or SPH+EoS in its (36, or 37+38, or 56), or in 
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Figure 13. Particle density, thermal energy per unit mass (in 
10^ units) and velocity for SPH and SPH+EoS modelling. Time 
T = 100. SPH+EoS [div(v)] refers to SPH+EoS modelling where 
eq. (60) is applied. 

its (60) formalism (hereinafter SPH+EoS [div(v)]) compare 
with each other and/or compare with the analytical solu- 
tion, taking into account that either the artificial viscosity 
in SPH or the physical dissipation in the SPH+EoS formal- 
ism (eqs. 36, or 37+38, or 56) activate whenever two par- 
ticles approach each other. At the same time we also check 
whether the SPH+EoS[div(v)] in its form (eq. 60) correctly 
behaves because its physical dissipation depends on V • t). 
Whichever is the final computational time, both v and the 
thermal energy per unit mass e for each particle should be 
constant. Any deviation from their initial values is an eval- 
uation of the inadequacy of the numerical approach to the 
true solution. Of course, we pay attention that, according to 
eqs. (42) and (43), the physical dissipation in the SPH+EoS 
formalism, according to eqs. (36, or 37+38, or 56, not 60), 
also depends on mutual particle separation and, to this pur- 
pose, the shorter is the vertical separation among particle 
rows, the smaller is the physical dissipation according to the 
Riemann problem solution, up to its natural limit towards 
the ID Riemann problem solution. 

Fig. 13 shows results of such simulations at time T — 
100 concerning particle density, thermal energy per unit 
mass and velocity (subsonic in such an example, with 7 = 
5/3 and e « 1.858-10^). Both SPH and, especially SPH+EoS 
(referred to eqs. 36, or 37+38, or 56, not 60), show statistical 
noise around v — 0, due to the fact that dissipation - artifi- 
cial, or physical - activates at the particle approaching and 
in particular among the flow particle in the pipe with the two 
tight boundaries, preventing the correct uniform slipping of 
the particle flow, since from the beginning of simulation. In- 
stead, results for the SPH+EoS case where the EoS takes 
into account of V • w (eq. 60), perfectly fit the analytical so- 
lution without any - or very small - vibrations in the velocity 
field, and consequently on e. Notice that such vibrations in 
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the velocity field are a bit larger when SPH+EoS is consid- 
ered according to eqs. (36, or 37+38, or 56, not 60). This 
is mainly due to the fact that the a and P counterpart are 
statistically larger than a = 1 and a = 1 + 2, as currently 
adopted in SPH for the solution of the Riemann problem. 
Density, as well as the thermal energy per unit mass are 
statistically not perturbed because of V • « even where 
w « both in SPH and in SPH+EoS approach. 



7 DISCUSSION AND CONCLUSIONS 

The comparison of ID Sod shock tube analytical results with 
SPIf-Riemann (SPH+EoS) ones, where a revision of the EoS 
within the Riemann problem is made, are quite successful. 
In particular, in our modelling, neither a parametrized ar- 
tificial viscosity, nor any dependence on spatial smoothing 
resolution length h, nor a sophisticated Godunov solver, nor 
additional computational time are involved. 

The simple EoS for inviscid ideal gases: p = (7 — l)pe 
cannot be strictly applied in shock problems because of the 
fact that, solving the Euler equations, not only instabili- 
ties and spurious heating come out, but also that this EoS 
derived by either the Charles, Volta Gay-Lussac and Boyle- 
Mariotte laws or by a partition function of statistical ther- 
modynamics, should be strictly applied only either in equi- 
librium configurations or to "reversible or quasi-static" evo- 
lutions of thermodynamic systems without any dissipation. 
This is a restriction that does not match with shock flow 
dynamics, when dissipation on the shock front cannot be ne- 
glected. Those techniques involving sophisticated Godunov- 
type schemes also introduce so me useful dissipation mech- 
anism in the numerical scheme l|Park fc Kwonll2003l ). right 
for obtaining correct solutions of Euler equations. 

To conclude, the general EoS here proposed, shows the 
correct behaviour, even in the presence of dissipative shocks 
in non viscous flows. A successful check of the reliability 
of the Riemann approach, here proposed on inviscid hydro- 
dynamics, is shown in App. A, where a study is shown on 
the coming out of spiral structures and of shocks in the ra- 
dial flow of accretion discs. The comparison, according to 
the different techniques, where artificial or physical dissi- 
pation is explicitly introduced, shows that the entire disc 
structure an d dyna mics, as wel l as the coming out of spi- 



rals (Sawadaetal 



ISawada fc Matsudal 



1987 : Spr uit et all 1 19871: iKaisiS 1 19891 : 



19921 : ISavonniie et al.lll994l ). are much 



better evidenced in those simulations where the EoS and 
the related dissipation are treated in their full physical sense 
according to the Riemann problem solution, where the dis- 
sipation is the most effective, while the EoS (eq. 60), where 
physical dissipation is correlated to the flow compression 
(V ■ 1;), shows the lower shear dissipation, being the particle 
tangential kinematics in the disk bulk more comparable to 
the Keplerian velocity fleld. 

T.l Concluding remarks 

• The need to introduce a numerical dissipation (implicit 
or explicit) is necessary to solve the hyperbolic Euler equa- 
tion system in non viscous flow dynamics. However, there is 
also a physical motivation because shock phenomena have 
to be considered as irreversible events. 



• Whichever is the adopted numerical technique for the 
non viscous hydrodynamics, a problem arises as far as 
the choice of some either explicit or implicit (adopting a 
Godunov-like technique) free parameters for dissipation is 
concerned. As an example, we recall the choice for a and /3 
artificial viscos ity coefficients in SPH, a*, (5 and Si in the 
formulation of lMorris fc MonaghanI ll 19971). th e calculation 
of pressure p* according to (|Monaghanlll997l ). as well as 
whet her either a the rmal energy diffusion term or an atten- 
uator iBalsaral l|l995l ) or both have to be used. 

• For these reasons, even if working in SPH, we propose 
an EoS for non viscous ideal gases that: 

introduce a physical dissipation both for the resolution 
of the Riemann problem - to solve shocks - and for the res- 
olution of shear flows; 

such a dissipation does not depend on arbitrary param- 
eters, as well as on spatial smoothing resolution length h. It 
is correlated to a and /3 through eqs. (42, 43, strictly valid 
for the Riemann problem), or through eqs. (62, 63 for a 
more general thematics) , according to local thermodynamic 
conditions; 

it is justified according to statistical thermodynamics 
calculations; 

• Analytic solutions for inviscid shocks and blast waves 
are highly idealized and subjected to Rankine-Hugoniot left- 
right {I, r) "jump conditions" limits {pi/pr ^ (7+l)/(7~l))- 
On theses conditions, analytical profiles are always "ruler- 
drawn". If shocks are irreversible thermodynamic events, a 
physical dissipation should smooth out every vertical linear 
profile of solutions; 

• In this work, ID profile of shocks and blast waves, where 
the physical dissipation is introduced in the EoS, better com- 
pare with the analytical solution and do not suffer from in- 
stabilities ("blimp") in the "flat" zones close to vertical dis- 
continuities in the physical parameters. Vertical descending 
profiles also better compare. Zones where rarefaction waves 
exist compare with SPH solutions. A larger discrepancy is 
evident only in the thermal energy profiles wherever verti- 
cal ascending profiles occur. However, we recall that thermal 
dissipation is not applied; 

• The equivalent formulations for the EoS of ideal gases 
(36, or 37+38, 56) appear successful for the treatment of 
the Riemann problem. Eq. (60) is also fair for non viscous 
shear flows. Nevertheless, even in not reflned form (36, or 
37+38 or 56), the physical dissipation, introduced in the 
EoS, better flts to the analytical sol ution of viscous flows if 
simply V Pi Cs h iMolteni et al.lll991^ describes the SPH dis- 
sipation (§6.3). The EoS in the form of eq. (60) introduces a 
physical dissipation only when a real local gas compression 
occurs. Hence, for accretion problems, if the velocity fleld 
is mainly Keplerian, either a real physical turbulent viscos- 
ity (|Flebbe et al.lll994l : lLanzafamell200^ . |2008| . |2009| ) in the 
Navier-Stokes approach, or an EoS in the Euler approach in 
the form of eqs. (36, or 37+38, or 56) should be considered; 

a and P (eqs. 42, 43) counterpart of the SPH+EoS in 
the form of eqs. (36, or 37+38, or 56) are generally larger 
than those values largely adopted in SPH techniques ~ 1. 
This lets us think that a = 1 and /? = 1+2 are a compromise, 
allowing description of both a Riemann problem well 
as a shear flow at the same time, with the price of including 
of some small shortcomings. The same a and values are 
mostly close to for those flows where no gas compression 
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occurs, being V ■ v ^ 0, when the SPH+EoS[div(v)] for 
of the Eos is taken into account in its more general (eq. 
60) form. However in some cases, wherever and whenever a 
high density increase occurs: p~^dp/dt = —V • i; 3> 0, both 
parameters could reach considerable values. 



APPENDIX A: THE APPEARANCE OF 
SPIRALS IN ACCRETION DISCS IN CLOSE 
BINARIES 

As an astrophysical application of EoS modification in fluid 
dynamics for a better identiflcation of shock profiles, we 
pay attention to the development of spiral patterns in ac- 
cretion discs in close binaries (CB), where collisions, shear 
fiows, turbulence and vorticity are concurrent throughout 
the structure and dynamics of the flows. Some remarkable 
highlights are here written on this theme to show whether, 
how and why the accuracy in the Riemann problem can af- 
fect the quality of results in non viscous flows. 

Results on comparison of 3D SPH simulations of accre- 
tion discs in CB are here shown as a further check with the 
aim of showing which method, including an explicit dissi- 
pation inviscid term in the EoS, gives results where clear 
spiral patterns and shock fronts in the radial flow are evi- 
dent. A rich s cientific literature exists on spiral patterns ap- 



pearance e.g. ijSawada e t al. 11198^1 ISpruit et al.l 



ll989l : ISawada fc Matsu da 1992; Savonniie et al. 



19871: Kaisii 
1 19941 ) 



Each 

simulation is stopped when a steady state configuration is 
obtained: i.e. when the number of particles within the grav- 
itational potential well of the primary compact star (e.g. 
a white dwarf or a neutron star) is statistically constant. 
The secondary star (a subgiant or a normal star) fills up 
entirely its Roche lobe, transferring its mass to the primary 
through the inner Lagrangian point LI. Disc's edges are free. 
iLanzafame et al.l (|2000l . 1200 il l clearly showed which geomet- 
ric and kinematic conditions favour the appearance of spi- 
ral structures in SPH. Therefore, in this section we refer to 
those initial (at LI point) and boundary conditions to high- 
light such structures. Particles move within the primary's 
Roche lobe towards the primary star at disc's inner edge 
(towards a central sphere whose radius equals 10~^), while 
they freely move outward from the disc's outer edge. Within 
the 10~^ radius they are considered in free fall, thus they 
are accreted. The smoothing resolution length of SPH parti- 
cles is ft = 5 ■ 10"'^, being 1 the non-dimensional separation 
rfi2 ~ 10® Km of the two stars. The injection of particles 
from LI is supersonic: Vinj = 130 Km , whilst the tem- 
perature of gas coming from the secondary star is T" = 10* K 
and 7 = 1.01. The compact primary is a IMq star, while 
the donor companion is a 0.5Mq star. 

The adoption of supersonic mass transfer conditions 
from LI is fully discussed in Lanzafame (2008, 2009), where 
disc instabilities, responsible for disc active phases of CB are 
discussed in the light of local thermodynamics. Whenever a 
relevant discrepancy exists in the mass density across the 
inner Lagrangian point LI between the two stellar Roche 
lobes, a supersonic mass transfer occurs as a consequence of 
the moment um flux conservation . The same result can also 
be obtained (|Lubow fc Shu|[l975l 'l by considering either the 
restricted problem of three bodies in terms of the Jacobi 
constant or the Bernoulli's theorem. 
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Figure Al. XY plots of 64 greytones density p isocontours of 
3D disc modelling in CBs. A linear scale of greytones is repre- 
sented on the left side, while its logarithmic scale is represented 
on the right. SPH refers to normal SPH results; SPH+EoS refers 
to SPH results where a reformulation of the EoS according to the 
Rieman problem solution is proposed; SPH+EoS [div(v)] refers to 
SPH results where a reformulation of the EoS in its more gen- 
eral physical sense is here proposed. When in statistically steady 
state, particles are ~ 79696 for SPH, ~ 72023 for SPH-l-EoS ans 
~ 112107 for SPH-(-EoS[div(v)]. 
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Figure A2. Radial distribution of deviation from the Keplerian 
kinematics, here represented as lj^—i^k, for the inner disc regions. 
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Figure A3. Radial distribution of the minimum and the max- 
imum values of a and /3 for each particle computed according 
to eqs. (42) and (43), respectively. Particles belonging to six 
Ar = 0.1 radial shells from r = 0tor = 0.6 are represented. 
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Figure A4. Radial distribution of the minimum and the max- 
imum values of a and /3 for each particle computed according 
to eqs. (62) and (63), respectively. Particles belonging to six 
Ar = 0.1 radial shells from r = to r = 0.6 are represented. 



The coming out of spiral patterns, as well as disc ellip- 
ticity are strongly evident in the disc model "SPH-f-EoS" 
of Fig. Al, where the EoS is reformulated according to 
a physical dissipation right for Riemann problem solution. 
The "SPH" disc model shows a glimmer of spiral structures, 
sti mulated by super s onic i njecti on conditions as evidenced 
in iLanzafame et all (|2000l . 1200 ll ). Instead, the plot of Fig. 



Al, relative to the disc model "SPH-fEoS[div(v)]", taking 
into account the full reformulation of the ideal gas equation 
in its non viscous physical sense, does not show such spiral 
characteristics. A central huge gas concentration, at a ra- 
dial distance corresponding to the angular momentum con- 
servation of particles injected from Ll, is shown. Hence, the 
whole disc's structure and kinematics is dominated by a cen- 
tral toroidal Keplerian structure since a physical dissipation, 
ruled out by local compression of gas bubbles, is rarely acti- 
vated. Ifowever, when a neighbour deviates from the strict 
Keplerian kinematics, it could be very strong. This demon- 
strates that the introduced physical dissipation within the 
EoS, suitable for the solution of the Riemann problem (eqs. 
36, or 374-38, or 56), affects the particle kinematics, pre- 
venting a too strong circularization of gas orbits in strictly 
Keplerian orbits. Such an "Eulerian" physical dissipation, as 
well as any other physical viscous dissipation, coming from a 
Navier-Stokes viscous, is absolutely necessary if the appear- 
ance of spiral structures are the desired dominant charac- 
teristics of disc structure. The stronger the dissipation, the 
higher the contrast, within some limitations. The prevention 
of any spiral has been noted in .Lanz afamc (2010), modelling 
the same accretion discs, where an arbitrary attenuation of 
artificia l viscosity, mai nly on a, has been applied accordin g 
both to lBalsaral l|l995h and to lMorris fc MonaghanI (|l997l ). 

The connection of disc ellipticity with sp iral patterns is 
accurately discussed inlBis ikalo et al.l (|l998l ). where a third 
spiral pattern can also develop in some models. Typically 
the two main spirals are on opposite sides in the disc struc- 
ture. The first one is directly connected to the incoming 
flow stream from Ll, while the second one comes from the 
more elongated disc outer edge on the opposite side to the 
injected infiow. The newest third one comes from the more 
elongated disc outer edge close to the injected stream. Other 
previous high c ompressibility (low 7) SPH non viscous disc 
models in C B (iMolteni et al.lll99ll: ILanzafame et al]|l992l 



1993al . Il99bl . |2000| . I2OOI 
20031 ) did not produce 



Yukawa et al. 1997 : Lanzafamg 



a such elliptical disc geometry. 
Flow perturbations, as well as tidal torques, are respon- 
sible for spiral appearance in disc structures. A very 
accurate discussion on tidal torque, its role and limits 



can b e found in iPaczvnskvl 1 19771): Papaloizou fc Pringld 
[ziiang fc ChenI (Il99^: llchikawa fc Osakil (|l992l . 



Lanzafamd (|2003l ). lii TLanzafamd l|2003l ) we paid at 



tention to the necessity to develop a Riemann-SPH code able 
to verify if weak shock fronts can develop within the disc 
bulk. After some time, we are now respecting that promise 
by considering a physical solution to the problem without 
any specific mathematical Riemann solver technique. 

Fig. A2 shows the radial distribution of the deviation 
of angular velocity component ujz to the Keplerian angu- 
lar velocity: — ojk in the inner disc regions. Of course, 
the "SPH-|-EoS[div(v)]" model shows the smaller deviation, 
both in its statistical sense (thickness of the distribution), 
and in the innermost radial regions, where the other two disc 
models progressively deviate, showing a lower angular mo- 
mentum component and, therefore a higher accretion rate 
onto the compact primary star. Together with Fig. Al, this 
further result shows, without any doubt, that the EoS per- 
fect gas reformulation (eq. 60) is really general because the 
"Eulerian" physical dissipation introduced is effective both 
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in the shear flow fluid dynamics and in the Ricmann problem 
flows at the same time. 

Fig. A3 shows the radial distribution of the minimum 
and the maximum values of a and /3, respectively, according 
to eqs. (42, 43), if the EoS in its (36, or 37-38, or 56) forms is 
considered. Those extreme values are computed comparing 
a and /3 for each neighbour for each particle. Both mini- 
mum and maximum a and /3 values converge towards a sin- 
gle value « 2 and « 1, respectively at LI due to the injector 
positioning. Sparc higher a and Rvalues for 0.5 ^ r ^ 0.6 
refer to other disc regions at its outer edge. In the disc bulk, 
the two distinct domains, on average, enlarge being in con- 
tact in the boundary of ~ 3.25 and ~ 2.75, respectively, 
since high temperature radial gradients characterize the in- 
ner disc zones, so that both "much colder" and "much hot- 
ter" neighbour companions affect the physical dissipation. 
All other intermediate a and /3 values, relative to each par- 
ticle (not represented for the sake of simplicity), span within 
the two separated domains for each panel of Fig. A3. Such 
plots show that a and /3 artificial viscosity counterpart of 
physical dissipation are larger, on the average, than the tra- 
ditional a = 1 and ^0 = 1 2, in so far as eqs. (36, or 37-1-38, 
or 56) arc taken into account. However, because any dissi- 
pation is physically prevented when rtj ■ Vij ^ 0, the a and 
P counterparts are not always effective. 

Fig. A4 shows, as in the previous figure, the radial dis- 
tribution of the minimum and the maximum values of a 
and /3, according to eqs. (62, 63), if the EoS in its (60) form 
is attributed to each particle. Even in this case, this deci- 
sion is taken to avoid too many points in the plots, repre- 
senting each neighbour particle. The scale in Fig. A4 is the 
same as for Fig. A3 for the sake of simplicity. Both a « 
and /3 ~ in the disc's outer regions, free of any gas com- 
pression, where local kinematics is mainly Keplerian at the 
disc's outer edge. This means that V • w is locally negli- 
gible. On the other external components, as the injected 
particle stream, no local compression occurs. For radial dis- 
tance r ^ 0.2, where the majority of disc's particles relies 
in the inner toroidal ring, both a and /S increase since the 
kinematics deviates from a perfect Keplerian behaviour and 
gas compression becomes effective there. This means that a 
non viscous hydrodynamics of a perfect gas, without a shear 
physical dissipation, also effective without any gas compres- 
sion, cannot produce any graduality in the radial transport. 
This implies that in an ideal gas, either the turbulent phys- 
ical dissipation in a Navier-Stokes approach, or an inviscid 
Eulerian fluid dynamics, where the Riemann problem is cor- 
rectly solved, must be considered to produce both a grad- 
uality in the radial transport, and spiral shaped structures 
and a statistically correct Keplerian tangential kinematics 
in a disc. 
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